#ifndef CSUPBR_GEOMETRY_H
#define CSUPBR_GEOMETRY_H

#include <exception>
#include <cassert>
#include <limits>
#include <cmath>
#include <algorithm>
#include <ostream>

namespace csupbr
{
	using Pixel = unsigned char;

	using Float = float;
	constexpr Float Pi = 3.1415926535f;
	constexpr Float Pi2 = Pi * 2.f;
	constexpr Float Pi4 = Pi * 4.f;
	constexpr Float InvPi = 1.0f / Pi;
	constexpr Float InvPi2 = 1.0f / (2.0f * Pi);
	constexpr Float InvPi4 = 1.0f / (4.0f * Pi);
	constexpr Float PiOver2 = Pi / 2.0f;
	constexpr Float PiOver4 = Pi / 4.0f;
	constexpr Float InvGamma = 1.0f / 2.2f;
    constexpr Float Inv256 = 1.0f / 255.9f;
	constexpr float TO_RGB = 255.9999f;
	constexpr float TO_RGB_INV = 1.f / 255.9999f;
	constexpr Float MaxFloat = std::numeric_limits<Float>::max();
	constexpr Float MinFloat = std::numeric_limits<Float>::min();
	constexpr Float MachineEpsilon = std::numeric_limits<Float>::epsilon();
	const Float RayBias = 0.05f;

	constexpr Float gamma(int n) {
		return n * MachineEpsilon / (1 - n * MachineEpsilon);
	}

	constexpr inline Float radiansToDegrees(Float rad) {
		return rad * InvPi * 180.0f;
	}

	constexpr inline Float degreesToRadians(Float deg) {
		return deg / 180.0f * Pi;
	}

	template<class T, size_t N>
	struct VecBase {
		static constexpr size_t _N = N;
		using _T = T;

		T _v[_N];

		VecBase() {
			for (auto& i : _v) {
				i = T();
				assert(!std::isnan(i));
			}
		}
	};

	template<class T>
	struct VecBase<T, 2> {
		static constexpr size_t _N = 2;
		using _T      = T;

		union {
			T _v[_N];
			struct {
				T x, y;
			};
			struct {
				T w, h;
			};
			struct {
				T u, v;
			};
		};

		VecBase() : x(T()), y(T()) {
		}

		VecBase(const T& x, const T& y) : x(x), y(y) {
			assert(!std::isnan(x) && !std::isnan(y));
		}

	};

	template<class T>
	struct VecBase<T, 3> {
		static constexpr size_t _N = 3;
		using _T      = T;

		union {
			T _v[_N];
			struct {
				T x, y, z;
			};
			struct {
				T r, g, b;
			};
			struct {
				// vertex, texCoords, normal
				T v, t, n;
			};
		};

		VecBase() : x(T()), y(T()), z(T()) {
		}

		VecBase(const T& x, const T& y, const T& z) : x(x), y(y), z(z) {
			assert(!std::isnan(x) && !std::isnan(y) && !std::isnan(z));
		}

	};

	template<class T>
	struct VecBase<T, 4> {
		static constexpr size_t _N = 4;
		using _T      = T;

		union {
			T _v[_N];
			struct {
				T x, y, z, w;
			};
			struct {
				T r, g, b, a;
			};
			struct {
				T i, j, k, l;
			};
		};

		VecBase() : x(T()), y(T()), z(T()), w(T()) {
		}

		VecBase(const T& x, const T& y, const T& z, const T& w)
			: x(x), y(y), z(z), w(w) {
			assert(!std::isnan(x) && !std::isnan(y) &&
				   !std::isnan(z) && !std::isnan(w));
		}

	};

	template<class T, size_t N>
	struct Vec : VecBase<T, N> {
		using VecBase<T, N>::VecBase;
		static constexpr size_t _N = VecBase<T, N>::_N;

		Vec& operator+=(const Vec& vec) {
			for (size_t i = 0; i < _N; i++) {
				this->_v[i] += vec._v[i];
			}
			return *this;
		}

		Vec& operator-=(const Vec& vec) {
			for (size_t i = 0; i < _N; i++) {
				this->_v[i] -= vec._v[i];
			}
			return *this;
		}

		Vec& operator*=(const Vec& vec) {
			for (size_t i = 0; i < _N; i++) {
				this->_v[i] *= vec._v[i];
			}
			return *this;
		}

		Vec& operator/=(const Vec& vec) {
			for (size_t i = 0; i < _N; i++) {
				this->_v[i] /= vec._v[i];
			}
			return *this;
		}

		Vec& operator*=(const T& vec) {
			for (size_t i = 0; i < _N; i++) {
				this->_v[i] *= vec;
			}
			return *this;
		}

		Vec& operator/=(const T& vec) {
			for (size_t i = 0; i < _N; i++) {
				this->_v[i] /= vec;
			}
			return *this;
		}

		Vec operator+(const Vec& vec) const {
			Vec tmp = *this;
			tmp += vec;
			return tmp;
		}

		Vec operator-(const Vec& vec) const {
			Vec tmp = *this;
			tmp -= vec;
			return tmp;
		}

		Vec operator*(const Vec& vec) const {
			Vec tmp = *this;
			tmp *= vec;
			return tmp;
		}

		Vec operator/(const Vec& vec) const {
			Vec tmp = *this;
			tmp /= vec;
			return tmp;
		}

		Vec operator*(const T& val) const {
			Vec tmp = *this;
			tmp *= val;
			return tmp;
		}

		Vec operator/(const T& val) const {
			Vec tmp = *this;
			tmp /= val;
			return tmp;
		}

		Vec abs() const {
		    Vec ret;
            for (int i = 0; i < N; i++) {
                ret[i] = std::abs(this->_v[i]);
            }
            return ret;
		}

		T dot(const Vec& vec) const {
			T sum = this->_v[0] * vec._v[0];
			for (size_t _i = 1; _i < N; _i++) {
				sum += this->_v[_i] * vec._v[_i];
			}
			return sum;
		}

		T dot() const {
			T sum = this->_v[0] * this->_v[0];
			for (size_t i = 1; i < N; i++) {
				sum += this->_v[i] * this->_v[i];
			}
			return sum;
		}

		T absDot(const Vec& v) const {
			return std::abs(dot(v));
		}

        Vec faceForward(const Vec &v2) {
            return (this->dot(v2) < 0.f) ? -*this : *this;
        }

		T lengthSquared() const {
			return dot(*this);
		}

		T length() const {
			return std::sqrt(lengthSquared());
		}

		void normalize() {
			float len = length();
			(*this) /= (len > std::numeric_limits<float>::epsilon()) ? len : MinFloat;
		}

		Vec<T, N> normalized() const {
			auto val = *this;
			val.normalize();
			return val;
		}

		Vec<T, N> operator-() const {
			auto tmp = *this;
			for (auto i = 0; i < N; i++) {
				tmp._v[i] = -tmp._v[i];
			}
			return tmp;
		}

		const T& operator[](size_t index) const {
			return this->_v[index];
		}

		T& operator[](size_t index) {
			return this->_v[index];
		}

		T max() const {
			T v = this->_v[0];
			for (int index = 1; index < N; index++) {
				v = std::max(v, this->_v[index]);
			}
			return v;
		}

		T min() const {
			T v = this->_v[0];
			for (int index = 1; index < N; index++) {
				v = std::min(v, this->_v[index]);
			}
			return v;
		}

		friend Vec<T, N> operator*(T k, const Vec<T, N>& rhs) {
			auto t = rhs;
			t *= k;
			return t;
		}

		static Vec lerp(const Vec& a, const Vec& b, const Vec& x)
		{
			Vec vx = b - a;
			return vx * x;
		}

		T distance(const Vec& vec) const {
			T sum = (this->_v[0] - vec._v[0]) * (this->_v[0] - vec._v[0]);
			for (size_t index = 1; index < N; index++) {
				sum += (this->_v[index] - vec._v[index]) * (this->_v[index] - vec._v[index]);
			}
			return sum;
		}

		static T distance(const Vec& vec1, const Vec& vec2) {
			T sum = (vec1._v[0] - vec2._v[0]) * (vec1._v[0] - vec2._v[0]);
			for (size_t i = 1; i < N; i++) {
				sum += (vec1._v[i] - vec2._v[i]) * (vec1._v[i] - vec2._v[i]);
			}
			return sum;
		}

		static Vec min(const Vec& a, const Vec& b) {
			Vec re;
			for (size_t i = 0; i < N; i++) {
				re[i] = std::min(a[i], b[i]);
			}
			return re;
		}

		static Vec max(const Vec& a, const Vec& b) {
			Vec re;
			for (size_t i = 0; i < N; i++) {
				re[i] = std::max(a[i], b[i]);
			}
			return re;
		}

	};

	template<typename T, size_t N>
	inline Vec<T, N> normalize(const Vec<T, N>& v) {
		return v.normalized();
	}

	template<typename T, size_t N>
	inline T dot(const Vec<T, N>& v) {
		return v.dot();
	}

	template<typename T, size_t N>
	inline T dot(const Vec<T, N>& a, const Vec<T, N>& b) {
		return a.dot(b);
	}

	struct Normal : Vec<float, 3> {
	public:
		Normal(const Vec<float, 3>& v) : Vec(v.x, v.y, v.z) {}

		Normal(float x, float y, float z) : Vec(x, y, z) {}

		Normal(float v = 0) : Vec(v, v, v) {}

	};

	using Point2i = Vec<int, 2>;
	using Point3i = Vec<int, 3>;

	using Vec2i = Vec<int, 2>;
	using Vec3i = Vec<int, 3>;

	using Point2f = Vec<float, 2>;
	using Point3f = Vec<float, 3>;

	using Vertex = Vec<int, 3>;

	struct Vec4f : Vec<float, 4> {

		Vec4f(const Vec<float, 4>& v) : Vec(v.x, v.y, v.z, v.w) {}

		Vec4f(float x, float y, float z, float w) : Vec(x, y, z, w) {}

		Vec4f(float v = 0) : Vec(v, v, v, v) {}

	};

	struct Vec3f : Vec<float, 3> {

		Vec3f(const Vec<float, 3>& v) : Vec(v.x, v.y, v.z) {}

		Vec3f(float x, float y, float z) : Vec(x, y, z) {}

		Vec3f(float v = 0) : Vec(v, v, v) {}

		Vec3f cross(const Vec3f& v) const {
			return Vec3f(
				y * v.z - z * v.y,
				z * v.x - x * v.z,
				x * v.y - y * v.x
			);
		}

		operator Vec4f() {
			return Vec4f(x, y, z, 0);
		}

	};

	inline Vec3f cross(const Vec3f& a, const Vec3f& b) {
		return a.cross(b);
	}

	struct Vec2f : Vec<float, 2> {
		Vec2f(const Vec<float, 2>& v) : Vec(v.x, v.y) {}

		Vec2f(float x, float y) : Vec(x, y) {}

		Vec2f(float v = 0) : Vec(v, v) {}

		operator Vec3f() {
			return Vec3f(x, y, 0);
		}

	};

	inline Float clamp(Float val, Float min_val, Float max_val) {
		return std::clamp(val, min_val, max_val);
	}

	inline Vec<float, 2> clamp(const Vec<float, 2>& p, const Vec<float, 2>& minp, const Vec<float, 2>& maxp) {
		return Vec<float, 2>{
			clamp(p.x, minp.x, maxp.x),
			clamp(p.y, minp.y, maxp.y),
		};
	}

	inline Vec<float, 3> clamp(const Vec<float, 3>& p, const Vec<float, 3>& minp, const Vec<float, 3>& maxp) {
		return Vec<float, 3>{
			clamp(p.x, minp.x, maxp.x),
			clamp(p.y, minp.y, maxp.y),
			clamp(p.z, minp.z, maxp.z),
		};
	}

	struct Mat4 {
	public:

		Mat4() = default;

		explicit Mat4(float val) {
			for (size_t i = 0; i < 4; i++) {
				for (size_t j = 0; j < 4; j++) {
					_rows[i][j] = (i == j) ? val : .0f;
				}
			}
		}

		explicit Mat4(float data[4][4]) {
			for (size_t i = 0; i < 4; i++) {
				for (size_t j = 0; j < 4; j++) {
					_rows[i][j] = data[i][j];
				}
			}
		}

		static Mat4 identity() {
			return Mat4(1);
		}

		static Mat4 translate(const Vec3f& v) {
			float t[4][4] = {
				{1, 0, 0, v.x},
				{0, 1, 0, v.y},
				{0, 0, 1, v.z},
				{0, 0, 0, 1},
			};
			return Mat4(t);
		}

		static Mat4 scale(const Vec3f& v) {
			float s[4][4] = {
				{v.x, 0,   0,   0},
				{0,   v.y, 0,   0},
				{0,   0,   v.z, 0},
				{0,   0,   0,   1},
			};
			return Mat4(s);
		}

		static Mat4 rotate(const Vec3f& axis, float angle) {
			const Float s = std::sin(angle);
			const Float c = std::cos(angle);
			const Float oc = Float(1.0) - c;
			float r[4][4] = {
				{
					oc * axis.x * axis.x + c,
					oc * axis.x * axis.y - axis.z * s,
					oc * axis.z * axis.x + axis.y * s,
					0
				},
				{
					oc * axis.x * axis.y + axis.z * s,
					oc * axis.y * axis.y + c,
					oc * axis.y * axis.z - axis.x * s,
					0
				},
				{
					oc * axis.z * axis.x - axis.y * s,
					oc * axis.y * axis.z + axis.x * s,
					oc * axis.z * axis.z + c,
					0
				},
				{0, 0, 0, 1}
			};
			return Mat4(r);
		}

		static Mat4 lookAt(const Vec3f& from, const Vec3f& to) {
			Vec3f up(0, 1, 0);
			Vec3f d = to - from;
			d.normalize();
			Vec3f xAxis = d.cross(up).normalized();
			Vec3f yAxis = xAxis.cross(d).normalized();
			float m[4][4] = {
					{xAxis.x, yAxis.x, d.x, from.x},
					{xAxis.y, yAxis.y, d.y, from.y},
					{xAxis.z, yAxis.z, d.z, from.z},
					{0,       0,       0,   1}
			};
			return Mat4(m);
		}

		static Mat4 perspective(float fov, float near, float far) {
			float invTanAngle = 1.0f / tanf(degreesToRadians(fov) * 0.5f);
			float p[4][4] = {
				{invTanAngle, 0, 0, 0},
				{0, invTanAngle, 0, 0},
				{0, 0, far / (far - near), -far * near / (far - near) },
				{0, 0, 1, 0}
			};
			return Mat4(p);
		}

		// static Mat4 orthographic(float near, float far) {
		// 	return Mat4();
		// }

		inline Vec<float, 4> operator[](size_t i)const {
			return _rows[i];
		}

		inline Vec<float, 4> row(size_t i)const {
			return _rows[i];
		}

		Vec<float, 4> column(size_t i) const {
			return Vec<float, 4>{_rows[0][i], _rows[1][i], _rows[2][i], _rows[3][i]};
		}

		Mat4 operator*(const Mat4& mat) const {
			Mat4 m = *this;
			for (size_t i = 0; i < 4; i++) {
				for (size_t j = 0; j < 4; j++) {
					m._rows[i][j] = _rows[i].dot(mat.column(j));
				}
			}
			return m;
		}

		Mat4& operator*=(const Mat4& mat) {
			auto m = *this * mat;
			for (int i = 0; i < 4; i++)
				_rows[i] = m._rows[i];
			return *this;
		}

		Vec<float, 4> operator*(const Vec<float, 4>& v) const {
			return Vec<float, 4>{
				_rows[0].dot(v), _rows[1].dot(v),
					_rows[2].dot(v), _rows[3].dot(v)
			};
		}

		Mat4 inverse(bool* suc = nullptr) const {
			auto m = reinterpret_cast<const float*>(_rows);
			float inv[16], det;
			int i;

			inv[0] = m[5] * m[10] * m[15] -
				m[5] * m[11] * m[14] -
				m[9] * m[6] * m[15] +
				m[9] * m[7] * m[14] +
				m[13] * m[6] * m[11] -
				m[13] * m[7] * m[10];

			inv[4] = -m[4] * m[10] * m[15] +
				m[4] * m[11] * m[14] +
				m[8] * m[6] * m[15] -
				m[8] * m[7] * m[14] -
				m[12] * m[6] * m[11] +
				m[12] * m[7] * m[10];

			inv[8] = m[4] * m[9] * m[15] -
				m[4] * m[11] * m[13] -
				m[8] * m[5] * m[15] +
				m[8] * m[7] * m[13] +
				m[12] * m[5] * m[11] -
				m[12] * m[7] * m[9];

			inv[12] = -m[4] * m[9] * m[14] +
				m[4] * m[10] * m[13] +
				m[8] * m[5] * m[14] -
				m[8] * m[6] * m[13] -
				m[12] * m[5] * m[10] +
				m[12] * m[6] * m[9];

			inv[1] = -m[1] * m[10] * m[15] +
				m[1] * m[11] * m[14] +
				m[9] * m[2] * m[15] -
				m[9] * m[3] * m[14] -
				m[13] * m[2] * m[11] +
				m[13] * m[3] * m[10];

			inv[5] = m[0] * m[10] * m[15] -
				m[0] * m[11] * m[14] -
				m[8] * m[2] * m[15] +
				m[8] * m[3] * m[14] +
				m[12] * m[2] * m[11] -
				m[12] * m[3] * m[10];

			inv[9] = -m[0] * m[9] * m[15] +
				m[0] * m[11] * m[13] +
				m[8] * m[1] * m[15] -
				m[8] * m[3] * m[13] -
				m[12] * m[1] * m[11] +
				m[12] * m[3] * m[9];

			inv[13] = m[0] * m[9] * m[14] -
				m[0] * m[10] * m[13] -
				m[8] * m[1] * m[14] +
				m[8] * m[2] * m[13] +
				m[12] * m[1] * m[10] -
				m[12] * m[2] * m[9];

			inv[2] = m[1] * m[6] * m[15] -
				m[1] * m[7] * m[14] -
				m[5] * m[2] * m[15] +
				m[5] * m[3] * m[14] +
				m[13] * m[2] * m[7] -
				m[13] * m[3] * m[6];

			inv[6] = -m[0] * m[6] * m[15] +
				m[0] * m[7] * m[14] +
				m[4] * m[2] * m[15] -
				m[4] * m[3] * m[14] -
				m[12] * m[2] * m[7] +
				m[12] * m[3] * m[6];

			inv[10] = m[0] * m[5] * m[15] -
				m[0] * m[7] * m[13] -
				m[4] * m[1] * m[15] +
				m[4] * m[3] * m[13] +
				m[12] * m[1] * m[7] -
				m[12] * m[3] * m[5];

			inv[14] = -m[0] * m[5] * m[14] +
				m[0] * m[6] * m[13] +
				m[4] * m[1] * m[14] -
				m[4] * m[2] * m[13] -
				m[12] * m[1] * m[6] +
				m[12] * m[2] * m[5];

			inv[3] = -m[1] * m[6] * m[11] +
				m[1] * m[7] * m[10] +
				m[5] * m[2] * m[11] -
				m[5] * m[3] * m[10] -
				m[9] * m[2] * m[7] +
				m[9] * m[3] * m[6];

			inv[7] = m[0] * m[6] * m[11] -
				m[0] * m[7] * m[10] -
				m[4] * m[2] * m[11] +
				m[4] * m[3] * m[10] +
				m[8] * m[2] * m[7] -
				m[8] * m[3] * m[6];

			inv[11] = -m[0] * m[5] * m[11] +
				m[0] * m[7] * m[9] +
				m[4] * m[1] * m[11] -
				m[4] * m[3] * m[9] -
				m[8] * m[1] * m[7] +
				m[8] * m[3] * m[5];

			inv[15] = m[0] * m[5] * m[10] -
				m[0] * m[6] * m[9] -
				m[4] * m[1] * m[10] +
				m[4] * m[2] * m[9] +
				m[8] * m[1] * m[6] -
				m[8] * m[2] * m[5];

			det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];

			if (det == 0) {
				if (suc) {
					*suc = false;
				}
			}

			det = 1.0f / det;

			Mat4 out;
			auto invOut = reinterpret_cast<float*>(out._rows);
			for (i = 0; i < 16; i++)
				invOut[i] = inv[i] * det;
			if (suc) {
				*suc = true;
			}
			return out;
		}
		
		Vec<float, 4> _rows[4];
		static_assert(sizeof(_rows) == sizeof(float) * 16, "Matrix4 must have packed 16 floats");
	};
	
	struct RGB : Vec<float, 3> {
		RGB()
			:Vec(0, 0, 0)
		{}

		RGB(float r, float g, float b)
			:Vec(r, g, b)
		{}

		RGB(unsigned int rgb) {
#define TO_RGB 255.9999f
			constexpr float INV_TO_RGB = 1.0f / TO_RGB;
			r = ((rgb >> 16 & 0xff) * INV_TO_RGB);
			g = ((rgb >> 8 & 0xff) * INV_TO_RGB);
			b = ((rgb & 0xff) * INV_TO_RGB);
#undef TO_RGB

		}

	};

	struct RGBA : Vec<float, 4> {
		RGBA()
			:Vec(0, 0, 0, 1.0f)
		{}

		RGBA(float r, float g, float b, float a)
			:Vec(r, g, b, a)
		{}

		RGBA(unsigned int rgba) {
#define TO_RGB 255.9999f
			constexpr float INV_TO_RGB = 1.0f / TO_RGB;
			r = ((rgba >> 24 & 0xff) * INV_TO_RGB);
			g = ((rgba >> 16 & 0xff) * INV_TO_RGB);
			b = ((rgba >> 8 & 0xff) * INV_TO_RGB);
			a = ((rgba & 0xff) * INV_TO_RGB);
#undef TO_RGB

		}

	};

    struct Spectrum : public Vec3f {
        Spectrum(const Vec<float, 3> &v) : Vec3f(v.x, v.y, v.z) {}

        Spectrum(float x, float y, float z) : Vec3f(x, y, z) {}

        Spectrum(float v = 0) : Vec3f(v, v, v) {}

        inline bool is_black() {
            for (float i : _v)
                if (.0f != i)return false;
            return true;
        }
    };

	struct Ray {
		Point3f origin;
		Vec3f direct;
		mutable float mint, maxt;
		float time;

		Ray(const Point3f& origin = { 0,0,0 }, const Vec3f& direct = { 0,0,1.f }, float time = 0.0f)
			:origin(origin), direct(direct), mint(0.001f), maxt(MaxFloat), time(time)
		{}

		inline Point3f to(float t) const {
			return origin + t * direct;
		}

	};

	template<typename T, size_t N>
	struct BBox {
		Vec<T, N> pMin, pMax;

		explicit BBox(const Vec<T, N>& p = { MaxFloat,MaxFloat,MaxFloat })
			:pMin(p), pMax(-p)
		{}

		explicit BBox(const Vec<T, N>& pmin, const Vec<T, N>& pmax)
			:pMin(Vec<T, N>::min(pmin, pmax)), pMax(Vec<T, N>::max(pmin, pmax))
		{}

		static BBox Union(const BBox& a, const BBox& b) {
			return BBox{ Vec<T,N>::min(a.pMin, b.pMin), Vec<T,N>::max(a.pMax, b.pMax) };
		}

		static BBox Union(const BBox& a, const Vec<T, N>& p) {
			return BBox{ Vec<T,N>::min(a.pMin, p), Vec<T,N>::max(a.pMax, p) };
		}

		// is 2 box overlaped
		bool overlaps(const BBox& b) const {
			bool x = (pMax.x >= b.pMin.x) && (pMin.x <= b.pMax.x);
			bool y = (pMax.y >= b.pMin.y) && (pMin.y <= b.pMax.y);
			bool z = (pMax.z >= b.pMin.z) && (pMin.z <= b.pMax.z);
			return (x && y && z);
		}

		// is point inside
		bool inside(const Vec<T, N>& p) const {
			return (
				p.x >= pMin.x && p.x <= pMax.x &&
				p.x >= pMin.y && p.y <= pMax.y &&
				p.z >= pMin.z && p.z <= pMax.z
				);
		}

		// expand by single value (cube)
		void expand(float delta) {
			pMin -= Vec<T, N>(delta, delta, delta);
			pMax += Vec<T, N>(delta, delta, delta);
		}

		float surfaceArea() const {
			auto d = pMax - pMin;
			return 2.f * (d.x * d.y + d.x * d.z + d.y + d.z);
		}

		float volume() const {
			auto d = pMax - pMin;
			return d.x * d.y * d.z;
		}

		// max extent axis index
		int maxExtent() const {
			auto diag = pMax - pMin;
			if (diag.x > diag.y && diag.x > diag.z) return 0;
			if (diag.y > diag.z) return 1;
			return 2;
		}

		const Vec<T, N>& operator[](int i) const { return i == 0 ? pMin : pMax; }
		Vec<T, N>& operator[](int i) { return i == 0 ? pMin : pMax; }

		Vec<T, N> lerp(const Vec<T, N>& x) const {
			return Vec<T, N>::lerp(pMin, pMax, x);
		}

		Vec<T, N> offset(const Vec<T, N>& p) const {
			return {
				(p.x - pMin.x) / (pMax.x - pMin.x),
				(p.y - pMin.y) / (pMax.y - pMin.y),
				(p.z - pMin.z) / (pMax.z - pMin.z)
			};
		}

		void boundingSphere(Vec<T, N>& center, float& radius) const {
			center = .5f * pMin + .5f * pMax;
			radius = Inside(center) ? center.distance(center) : .0f;
		}

		bool intersectP(const Ray& ray, float& thit0, float& thit1) const {
			float tmin = ray.mint, tmax = ray.maxt;
			for (int i = 0; i < N; i++) {
				float invRayDir = 1.0f / ray.direct[i];
				float tNear = (pMin[i] - ray.origin[i]) * invRayDir;
				float tFar = (pMax[i] - ray.origin[i]) * invRayDir;
				if (tNear > tFar)std::swap(tNear, tFar);
				tmin = tNear > tmin ? tNear : tmin;
				tmax = tFar < tmax ? tFar : tmax;
				if (tmin > tmax)return false;
			}
			thit0 = tmin;
			thit1 = tmax;
			return true;
		}

	};

	using BBox2f = BBox<float, 2>;
	using BBox3f = BBox<float, 3>;

	template<typename T, size_t N>
	inline std::ostream& operator<<(std::ostream& os, const Vec<T, N>& v) {
		if (N == 1) os << "vec1[ ";
		if (N == 2) os << "vec2[ ";
		if (N == 3) os << "vec3[ ";
		if (N == 4) os << "vec4[ ";
		for (int i = 0; i < N; i++) {
			os << v._v[i];
			if (i != N - 1)
				os << ", ";
		}
		os << " ]";
		return os;
	}

	inline std::ostream& operator<<(std::ostream& os, const Mat4& mat) {
		os << "mat4[" << mat._rows[0]
		   << "\n     " << mat._rows[1]
		   << "\n     " << mat._rows[2]
		   << "\n     " << mat._rows[3] << "]";
		return os;
	}

	inline std::ostream& operator<<(std::ostream& os, const Ray& ray) {
		os << "ray[" << ray.origin << " -> " << ray.direct << "]";
		return os;
	}

	template<typename T, size_t N>
	inline std::ostream& operator<<(std::ostream& os, const BBox<T, N>& bbox) {
		if (N == 2) os << "bbox2[" << bbox.pMin << " >-< " << bbox.pMax << "]";
		if (N == 3) os << "bbox3[" << bbox.pMin << " >-< " << bbox.pMax << "]";
		return os;
	}

}

#endif // !CSUPBR_GEOMETRY_H